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ABSTRACT 

We present an integrated photometric spectral energy distribution (SED) of the 
Magellanic-type galaxy NGC 4449 from the far-ultraviolet to the submillimetre, in- 
cluding new observations acquired by the Herschel Space Observatory. We include 
integrated UV photometry from the Swift Ultraviolet and Optical Telescope using a 
measurement technique which is appropriate for extended sources with coincidence 
loss. 

In this paper, we examine the available multiwavelength data to infer a range 
of ages, nietallicities and star formation rates for the underlying stellar populations, 
as well as the composition and the total mass of dust in NGC 4449. Our analysis 
of the global optical spectrum of NGC 4449 fitted using the spectral fitting code 
STARLIGHT suggests that the majority of stellar mass resides in old (> 1 Gyr old) 
and metal-poor {Z/Zq ^ 0.2) populations, with the first onset of star formation activ- 
ity deduced to have taken place at an early epoch, approximately 12 Gyr ago. A simple 
chemical evolution model, suitable for a galaxy continuously forming stars, suggests 
a ratio of carbon to silicate dust mass comparable to that of the Large Magellanic 
Cloud over the inferred time-scales. 

We present an iterative scheme, which allows us to build an in-depth and multi- 
component representation of NGC 4449 'bottom-up', taking advantage of the broad 
capabilities of the photoionization and radiative transfer code MOCASSIN (MOnte 
CArlo Simulations of Ionized Nebulae). We fit the observed SED, the global ionization 
structure and the emission line intensities, and infer a recent star formation rate of 
0.4 M0yr~^ and a total stellar mass of « 1 x 10^ Mq emitting with a bolometric 
luminosity of 5.7 x 10^ L©. Our fits yield a total dust mass of 2.9 ± 0.5 x 10^ Mq 
including 2 per cent attributed to polycyclic aromatic hydrocarbons. We deduce a 
dust to gas mass ratio of 1/190 within the modelled region. While we do not consider 
possible additional contributions from even colder dust, we note that including the 
extended H i envelope and the molecular gas is likely to bring the ratio down to as 
low as ~ 1/800. 
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1 INTRODUCTION 

The conditions in the interstellar medium (ISM) in a galaxy 
and the environment-dependent processes governing the for- 
mation of stars and dust can be studied on many physical 
scales. Considering an entire galaxy, the interplay between 
the underlying stellar emission, the degree of radiation re- 
processing taking place in the ISM and the thermal or non- 
thermal emission due to dynamical or evolutionary processes 
give rise to a unique pattern of emission and absorption fea- 
tures in the observed spectral energy distribution (SED). 

The spectroscopic measurements alone provide a wealth 
of information about the ionized (Hll) regions or photo- 
dissociation regions (PDRs) in galaxies (e.g. Guseva et al. 
2004 Vasta et al.||20l6 l. Studies constrained by broadband 
photometric measurements can allow a better understanding 
of the main components of the thermal continuum, namely, 
the stellar emission and the emission due to interstellar 
dust. Observed emission lines and diagnostic line ratios pro- 



vide important further constraints (Martmez-Galarza et al. 



20111. Synthetic single stellar populations (SSPs) can be 



used to characterise star-forming regions within galaxies 
(e.g. Groves et al. 2008 Ferreras et al. 20121, or to de- 



compose the observed SEDs of entire galaxies (e.g. lAmormI 



et al. 20121. The composition of dust and the masses of 



dust species can be inferred from observations by assuming 



measured laboratory properties of dust (e.g. Galliano et al 



2008a and references therein). Together, these components 



can provide an in-depth view of a galaxy (e.g. Cormier et al 
20T2I [Hermelo et al.|20"T3| . 

Numerical models reproducing all the available photo- 
metric and spectroscopic measurements across a wide range 
of wavelengths can provide the most complete picture of the 
integrated properties of a galaxy. We wish to examine a mul- 
ticomponent model treating the stellar content, the gaseous 
phase and the dust within a galaxy in a self-consistent way. 
Our initial goals include finding a robust numerical scheme, 
which could be applied to a sample of dwarf galaxies in or- 
der to study the details of their star formation histories and 
dust content. 

In this paper, we describe the first multicomponent 
observation-driven model of a galaxy generated with the 
photoionization and radiative transfer code MOCASSIN 
l Ercolano et al.||2003| |2005| [2008] ), which includes a simul- 
taneous and self-consistent treatment of the stellar compo- 
nent, the gaseous phase, its ionization structure and mul- 
tiple dust species. MOCASSIN is a fully three-dimensional 
code, whose recent applications include modelling the ex- 
treme bipolar planetary nebula NGC 6302 ( [Wright et al 
2011 1 and studying dust emission by supernova ejecta (Wes- 



son et al. 2010). Applying MOCASSIN to simulate entire 



galaxies is more resource-consuming since both the gas phase 
and the dust phase need to be treated self-consistently over 
relatively large physical scales. The small intrinsic size of 
dwarf galaxies make them ideal candidates for such mul- 
tiwavelength studies. Indeed, dwarf galaxies provide small 
self-contained but dynamic environments. As such, they are 
useful in developing our understanding of mechanisms un- 
derlying star formation, including chemical enrichment and 
feedback processes, especially for low-metallicity objects re- 
sembling those formed at earlier epochs in the evolution of 
the Universe. 



In a similar context 
model Mrk 996 



MOCASSIN was first used to 
In this paper we describe 



20091 



a model of the well-known starburst galaxy NGC 4449, 
which is a low-metallici ty [log(0/H) -I- 12 = 8.23, 1/3Zq; 
Engelbracht et al. 2008 actively star-forming (SF R ~ 0.5 

barred 



M 



©yr 



Hunter et al. 



1998 



Hill et al. 



19981 



Magellanic-type irregular galaxy (de Vaucouleurs et al. 1991 
Corwin et aL]|1994 l seen face-on. At a distance of 3.8 Mpc 



(Annibali et al. 20081 its Hi envelope extends to 12.9 kpc 



(11.6), equal to approximately three Holmberg radii ( Swa- 
ters et al.|2002 Hunter et al.|1 99"9') , and forms a pronounced 



system of streamers with a counter-rotating H I core, which 
may be indicative of past mergers ( [Hunter et al.|1999[[Theis 
[& Kohle"2001|). NGC 4449 is particularly suitable for mul- 
tiwavolcngth studies because of its wide range of existing 
photometric data extending to the far-infrared (FIR). The 
new observations acquired with the Spitzer Space Telescope 



(Werner et al. 



(Pilbratt et al. 



20041 and the Herschel Space Observatory 



2010 \ are essential in studying gas cooling 



processes and the properties of dust in the ISM of NGC 4449. 

In the following sections we discuss in detail a model of 
NGC 4449 based on multiple observational constraints and 
including a treatment of polycyclic aromatic hydrocarbons 
(PAHs). Thus, in Section |2] we describe the observations 
and data reduction techniques, and in Section[3|we describe 
the modelling method, as well as the assumptions and the 
convergence criteria adopted. The results and discussion are 
presented in Section [4[ and the conclusions follow in Sec- 
tionlsl In a companion paper, Karczewski et al. (2013 here- 



r|5L 

jrTlF 



after. Paper III we study the FIR cooling lines and discuss 



properties of PDRs inferred from spatially-resolved spectro- 
scopic observations of NGC 4449. 



2 PHOTOMETRIC OBSERVATIONS AND 
DATA REDUCTION 

In constructing the observed SED of NGC 4449 we have 
used data from previous studies (eg, Kcnnicutt 1992, En-[ 
gelbracht et al.|2008||Bendo et al. 2012b, , Hunter et al. 19861 



Bottner et al. 2003 1 and from systematic surveys [Sloan Digi- 
tal Sky Survey (SDSS),[York et al.|2000' Two Micron All Sky 



Surve y (2MASS), [Skrutskie et al.|20 06 W/gE, [^AWght et al 
2010 and Planck, Planck Collaboration||2011 . We present 
archival data fro m GALEX ([Martin et al.[[2005} UV wave- 
bands) and Swift ( Gehrels et al.|2004 UV and optical wave- 
bands), as well as new FIR photometry from Herschel (Pil- 
bratt et al.[[2010] ). A journal of observations is given in Ta- 
bleiU 

In Table [2| we list all global photometric measurements 
used to construct the observed SED for NGC 4449. Ex- 
isting submillimetre SCUBA (Submillimetre Common-User 
Bolometer Array) measurements at 450 fim and 850 ^im 
( Bottner et al. 2003 1 were omitted from the list of photo- 



metric measurements due to their incomplete field-of-view 
(FOV) and the availability of space borne Herschel and 
Planck observations offering a similar wavelength coverage. 
Five illustrative false-colour images of NGC 4449 over se- 
lected broadband ranges are presented in Fig. [l] 
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2.1 Ultraviolet observations with GALEX and 
Swift 

At UV wavelengths NGC 4449 has been observed by two re- 
cent space-borne missions, GALEX (in two UV bands) and 
the Swift Ultraviolet /Optical Telescope (UVOT; in three 
UV and three optical bands) . Standard aperture photometry 



using the in-flight calibration of Morrissey et al. ( 2007 1 was 



performed on tile 5228 acquired as part of GALEX Nearby 
Galaxy Survey and accessible via General Release 6. The 
overall uncertainties include photometric repeatability mea- 



surements by Morrissey et al. ( 2007 1 



Swift /UVOT ( [Rorning et al.|2005| ) offers four times the 



angular resolution of GALEX in three narrow UV bands. 
The telescope is equipped with a photon-counting detector 
and was originally designed to detect and observe gamma- 
ray bursts. As a consequence, standard data reduction pro- 
cedures cannot be applied straightforwardly to observations 
of extended sources. In the Appendix we discuss this prob- 
lem in detail and present a method of obtaining global pho- 
tometry, which is appropriate for extended sources. 



2.2 Optical observations with SDSS 

In addition to the optical fluxes derived from the 
Swift /VYOT observations, SDSS observations in five opti- 
cal bands {u, g, r, i and z) are available as part of Data 
Release 7. NGC 4449 is contained entirely within one field, 
thus minimizing deblending or sky subtraction issues which 
may arise in the case of extended sources. Standard aperture 
photometry was performed using elliptical apertures and the 
counts were converted into fiux densities as described by 



West et al. (20101. The uncertainties in the derived fluxes 



are dominated by uncertainties in sky determination and 
subtraction. 



2.3 Far-infrared observations with Herschel 



NGC 4449 has been observed by Herschel (Pilbratt et al 



2010 ) as part of the Dwarf Galaxy Survey (DGS; Madden 
et al.|2013 l with its two imaging photometers: the Photode- 



tector Array Camera and Spectrometer (PACS; !Poglitsch' 
et al.|2010| ) at 70 fim, 100 fim and 160 fim and the Spectral 



and Photometric Imaging REceiver (SPIRE; Griffin et al 



2010) at 250 ^im, 350 ^m and 500 fim. The Fufl Width at 



Half Maximum (FWHM) of the PSF is 5'.'2, 7'.'7, 12'.'0, 18'.'2, 
24'.'9, 36'.'3 in these bands, respectively. The details of all 
PACS and SPIRE data reduction steps, including error es- 
timation, can be found in Remy-Ruyer et al. (20131. 



The PACS observations were performed as four pairs of 
orthogonal scans covering an area of 24' x 24'. We used an 
adapted version of the standard script of v7.0 of the Herschel 
Interactive Processing Environment (HIPE; Ott|2010 1. The 
basic processing includes fiagging bad or saturated pixels, 
converting the signal into Jy /pixel and applying the flatfield 
correction. Additionally, we systematically masked column 
of all of the constituent 16 x 16 matrices in the PACS array 
to avoid electronic crosstalk and we performed second level 
deglitching. The resulting Level 1 products were converted 
into maps with the pixel size of 2, 2 and 4 arcsec for the 



three PACS bands and processed with Scanamorphos ( Rous- 
sel|2012 1 , which is particularly suitable for extended sources 



with low frequency noise. The integrated PACS fluxes of 
NGC 4449 are well represented by a 40 K blackbody and 



suitable colour corrections were applied accordingly ( Miiller 



et al. 2011b I, where no additional correction factors are 



required to correct for the extended nature of the source 
(Sauvagc 2011). The quoted uncertainties include the cali- 



bration uncertainties at 5 per cent in all bands ( Remy-Ruyer 
et al.|2013| [Miiller et al.|2011a[ ) 



The SPIRE observations were performed as two orthog- 
onal scans covering an area of 24' x 24'. The corresponding 
maps were reduced using a modified version of the SPIRE 
pipeline in HIPE. The steps up to Level 1 were identical 
to those in the original pipeline provided by the SPIRE In- 
strument Control Center (ICC). Additionally, residual base- 
line subtraction was performed by subtracting the median 
of the time-lines for each bolometer over the entire observa- 
tion ( [Pohlen et aL]|2010[ ). This was followed by an iterative 
process to completely remove residual signals that appear as 
stripes in the maps ( Bendo et al.|2010 |. The final map was 
constructed using Naive Mapper available in HIPE and cal- 
ibrated for an extended source. A modified blackbody fit of 
the form S^ oc B^{T) yields /3 « 2 globally for NGC 4449 



(Remy-Ruyer et al. 20131, which was used for colour cor- 



rection of the integrated fluxes (SPIRE Observers' Man- 
ual; Valtchanov|2011 1. The uncertainties include the revised 
overall calibration uncertainties at 7 per cent in all SPIRE 
bands ( jCriflan fc Lim|2011| ). 



2.4 2MASS, WISE and Planck catalogue data 

Global 2MASS photometric measurements for NGC 4449 
were taken from columns j _m_ext, h_m_ext and k_m_ext of the 
2MASS All-Sky Extended Source Catalog, and were subse- 
quently converted to flux densities using the zero points tab- 



ulated by Cohen et al. ( 2003 1. WISE photometry was taken 



from columns gmag and gerr of the WISE All-Sky Source 
Catalog. Colour correction was applied using interpolated 
correction factors suitable for a source emitting as v~^'^^ 
across the four WISE bands. Planck photometry was taken 
from column GAUFLUX(_ERR) of the Early Release Compact 
Source Catalogue. Similarly, the interpolated colour correc- 
tion factors used were suitable for a source emitting as v^'^^ 
across the Planck bands centred at 350, 550 and 850 yLxn. 



3 MODELLING METHOD 
3.1 Overview 

The numerical code MOCASSIN (MOnte CArlo Simula- 



tionS of Ionized Nebulae; Ercolano et al. 2003 2005 20081 



was originally intended as a tool to construct realistic 
models of photoionized nebulae. It allows arbitrary three- 
dimensional geometries, separate for gas and dust, multiple 
ionizing sources emitting with a given input spectrum, vari- 
able gas chemistry, multiple dust species and arbitrary dust 
grain size distributions. Given these input parameters, the 
code self-consistently solves the radiative transfer in the co- 
existing gas and dust phases and calculates the ionization 
degree, electron temperature and dust grain temperature at 
every grid cell, and the overall emergent SED of the gas and 
dust. 
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Figure 1. False-colour (R, G, B) broadband images of NGC 4449: (a) Swift/VYOT (2486 A, 2221 A, 1991 A), (b) SDSS (r, 
g, u), (c) Spitzer/IRAC (8 ^lm, 5.8 fim, 3.6 ^lm), (d) Herschel /PACS (160 ^m, 100 /xm, 70 ^lm) and (e) Herschel/ SPIKE 
(500 fim, 350 fim, 250 /xm). North is up, east is to the left. The images are centred at 12'>28™11!1, +44°05'37" (J2000). 
The bar in the top left image is 2' in length (2.2 kpc). 



Instrument 


Bandpass/ Aoff 


Time (s) 


Date observed 


Observation ID 




New observations 






Herschel/PACS 


70 /jm + 160 lira 


3497 


2011 May 16 


1342221125 




70 /im + 160 lira 


3497 


2011 May 16 


1342221126 




100 /im + 160 /im 


3497 


2011 May 16 


1342221127 




100 /.tm + 160 /im 


3497 


2011 May 16 


1342221128 


Herschel/SPIRE 


250 ^ira, 350 /xm, 500 /.tm 


1035 


2010 June 12 


1342198243 




Archival observations 






Swift/VVOT 


/2030A 


1544 


2007 Mar 27 


00035873001 




TOm2/223lA 


1091 


2007 Mar 27 


00035873001 




uvwl /2634A 


771 


2007 Mar 27 


00035873001 




«/350lA 


385 


2007 Mar 27 


00035873001 




6/4329A 


385 


2007 Mar 27 


00035873001 




i;/5402A 


385 


2007 Mar 27 


00035873001 


GALEX 


FUV/1539A 


835 


2006 Mar 15 


tile 5228 




NUV/2316A 


765 


2006 Mar 15 


tile 5228 


SDSS 


■!i/355lA 


54 


2003 Mar 24 


3813/1/41/237/241 




t//4686A 


54 


2003 Mar 24 


3813/1/41/237/245 




r/6165A 


54 


2003 Mar 24 


3813/1/41/237/237 




i/748lA 


54 


2003 Mar 24 


3813/1/41/237/239 




z/893lA 


54 


2003 Mar 24 


3813/1/41/237/243 



Table 1. List of observations. 



In the sections that follow we will describe how MO- 
CASSIN (v2.02.70) can be used to model large physical sys- 
tems, such as galaxies. Our model is constructed 'bottom- 
up', where as many parameters as possible are fixed a priori 
based on observations. The simulations are set up with the 
empirical elemental abundances and the observed radial gas 
distribution fixed. The input stellar spectrum corresponding 
to a mixture of stellar populations, scaled by the stellar lu- 
minosity Li,, and the dust to gas mass ratio (A'/dust/Mgas) 
are free parameters. Based on these input parameters, which 



form a theoretical model of a galaxy, MOCASSIN produces 
a predicted low-resolution SED and a full set of predicted 
emission line intensities. If these predictions agree with ob- 
servations, the original parameter set can be regarded as 
a true representation of the galaxy under the assumptions 
made. Otherwise, individual parameters can be adjusted and 
simulations repeated iteratively, as summarised schemati- 
cally in Fig. [2] 

In constructing our MOCASSIN model of NGC 4449 
we have made simplifying assumptions about the geometry, 
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Survey/instrument 






Aperture" 


Reference 


GALEX 


1539 A 


152 ± 12 mjy 


11^5 


(1) 




2316 A 


183 ± 9 mJy 


7'. 5 




Swift/VYOT 


1991 A 


175 ± 12 mJy 


8'.5 X 6' 


(1) 




2221 A 


189 ± 14 mJy 








2486 A 


202 ± 15 mJy 








3442 A 


249 ± 22 mJy 








4321 A 


463 it 41 mJy 








O'^ilU A 


Dii It ( mjy 






SDSS 


3551 A 


242 it 5 mJy 


6'.5 X 4'7 


(1) 




4686 A 


457 it 14 mJy 








6165 A 


573 it 40 mJy 








i 4oi A 


Dzz It ou mjy 








8931 A 


683 it 34 mJy 






2MASS 


1.24 /xm 


916 it 21 mJy 


8 


(2) 




1.66 /xm 


1070 it 30 mJy 








2.16 /iiii 


ooy It oi mjy 






WISE 


3.4 /im 


391 it 7 mJy 


4'.0 X 2'8 


(3) 




4.6 fivci 


251 it 4 mJy 










yo4 It 1 ( mjy 








22 fiui 


2860 it 60 mJy 






Spitzer /IRAC 


3.6 /im 


493 it 15 mJy 




(4) 




4.5 /im 


ol i it lU mJy 








5.8 /^m 


615 it 19 mJy 








8 /^m 


1420 it 40 mJy 






i-tlAo 


12 /^m 


z.i it U.z Jy 


o 






25 


4.7 it 0.5 Jy 








60 ^lu 


36 it 3 Jy 






opilZcT f iVilr o 




o.zy It u. io Jy 


y.o X D.D 


\p) 




70 fira 


43.8 it 4.4 Jy 








160 /im 


78.1 it 9.4 Jy 






rievscnel / r^AL^o 


70 /^m 


4y.o it z.o Jy 


of n 
O.O 


{') 




100 


75.9 it 3.8 Jy 








160 /xm 


79.5 it 4.0 Jy 






Herschel/SmK& 


250 Atm 


35.7 it 2.5 Jy 


8'.3 


(7) 




350 /xm 


16.2 it 1.1 Jy 








500 Atm 


5.6 it 0.4 Jy 






Planck 


350 /xm 


16.14 it 0.34 Jy 


6'.1 X 5^2^= 


(8) 




550 /xm 


4.97 it 0.24 Jy 


5'.8 X 4^8^= 






850 /im 


1.45 it 0.15 Jy 


6'.0 X 4^4^= 




IRAM 


1.2 mm 


260 it 40 mJy 


two fields 2'.3 


(9) 



" given as diameters for circular apertures or major X minor axes for elliptical 
apertures; major axis of the 'total' aperture; FWHM of a Gaussian fit 

Table 2. Summary of the derived photometric data for NGC 4449. All fluxes have 
been corrected for the foreground extinction using EjB - V) = 0.019 l |Schlegel| 
|et al.|1998t and the extinction law of |CardeUi et al.| | |1989| l with Rv = 3.1. Ref- 
erences: (1) this work; (2) 2MASS All-Sky Extended Source Catalog | |Skrutskie 
et al.||2006t; (3 ) WIS E Al l-Sky Source Catalog jWright et al.||2010t; (4) |En- 
gelbracht et al.| ||2008|l; ( 5) [Hunter et al] | |1986| l; (6) |Bendo et al.| ( 12012^1 1; (7) 
Remy-Ruyer et al.| l|2013ll; ( 8) Planck Early Release Compact Source Catalogue 
planck Collaboration|2011^ ; (9) |Bottner et al.|p003l l. 



gas density distribution, elemental abundances, dust com- 
position and dust grain size distribution. 

We have also made assumptions about the star forma- 
tion history. In Section [374l we discuss a general picture of the 
stellar populations in NGC 4449, obtained using the spec- 
tral decomposition code STARLIGHT. In Section [3.7. 1| we 
characterise the youngest stellar population based on obser- 
vational constraints. Finally, in Section [3.7.3| we fit the older 



populations taking into account the constraints from earlier 
sections. 

The fitting of the star formation history and dust is 
performed almost simultaneously for self-consistency. There- 
fore, the young stellar population, the episodes of star forma- 
tion and dust are never modelled on their own. The youngest 
stellar population is fitted together with a representative 
sample of older populations, and the final dust content is 
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trial stellar 




spectrum 




trial L. 




trial Mdu8t/Mgas 





geometry, n(r) 

i 



Hi, optical. Ha 



MOCASSIN 



J' J' >l 



(iterate) 



low-res L(Ha) line ( compar e) broadband photometry, 
^ SED QCH"} fluxes ^ ^ global spectra 

i 

(accept) 

Figure 2. Summary of the modelling method. The shaded 
regions group together the input variables (left) the outputs 
(middle) and the observational constrains (right). 



determined relatively early in the process to ensure self- 
consistency. 

The assumptions listed below form basis of the iterative 
scheme described in detail in SectionFS.?! 



3.2 Gas density distribution 

Spherical symmetry is assumed. Although MOCASSIN is a 
fully three-dimensional code, this capability is not exploited 



at this time due to technical limitations (see Section 3.6 1. 
However, the method presented in this paper can be very 
easily adapted for fully three-dimensional modelling and the 
first step towards such modelling is discussed in Section 1X4] 



We used the H I radial profile of Swaters et al. ( 2002 1 to 
approximate nH(r)/cm~^ as four second-order polynomials 
in the form ar^ + br + c for radii up to 3.3 kpc (3'.0), cor- 
responding to the extent of the galaxy in the FIR (Fig. [T]). 
Although NGC 4449 is highly irregular, for simplicity we 
treat the star formation regions cumulatively, and based on 
the optical image (Fig. [T]) we assume that all gas is ionized 
within a distance of ~ 0.4 kpc from a single ionizing source. 
The total observed H i mass for NGC 4 449 is 2-2.5 x 10^ M© 
( Swaters et al.|2002 Bajaja et al. 19941, which is distributed 
between a system of streamers of mass 9 x 10* M© and 
the central regio n with diameter 24-35 kp c and mass 1.1- 

The H I mass 



1.25 X 10^ Mq (Hunter et al 



1998 



19991 



within the radius of 3.3 kpc derived using our formulation is 
5.5 X 10* Mq, which is consistent with the estimated total 
of 1.1-1.25 X 10^ Mq for the inner system given the radial 
surface density profile of [Swaters et al.| ( |2002[ ). 



The 



of H2, based on L, 



CO (1-01 



= 8.4 X 



10*^ K kms'^ pc^ (corrected for D = 3.8 Mpc; |Bottner et al 
2003[) and assuming qco ~ 10-20 M© (Kkms 'pc^)"" (cf. 



Sandstrom et al.|[20T2l [Schruba et"ar]|20T2 



and references 

therein), can be estimated at ~ 0.8-1.7 x 10* Mq. However, 
we note that CO-to-H2 conversion factors for low metallic- 
ity objects are very uncertain, and therefore we include no 
estimates of H2 mass in our gas mass totals. 

We assumed that all gas exists in small clumps de- 
scribed by the filling factor e, which is the same at all radii. 
Considering a Stromgren sphere of radius 0.4 kpc with Lua 



given by Hunter et al. ( 1999 1 and assuming that the average 
electron density is equal to the central ric given by IMartm] 
( 19971, we derive an initial estimate of e ~ 0.003. This is re- 



formalism of Barlow (19871 we estimate that the mass of 
ionized gas is ~ 1.5 X iFMq. 



3.3 Elemental abundances 

Elemental abundances are assumed to be constant through- 
out the galaxy. The abundances were taken from |Vigroux] 
et al. (19871, except for carbon, neon and sulphur, where 
the Large MageUanic Cloud (LMC) ratios for C/O, Ne/0 



and S/0 were assumed ( Russell fc Dopita| |1992'). Although 
the sophisticated framework offered by MOCASSIN can also 
be used to determine elemental abundances ( jErcolano et al.| 



20101, in this paper MOCASSIN is not used in this context. 



3.4 Star formation history 

Star formation histories are difficult to constrain because of 
observational limitations (e.g. instrumental sensitivity), but 
also because of the more fundamental age-metallicity de- 
generacy (e.g. |Worthey]|1994[ [Ferreras et al.|[T999[ |Ferrer£is| 
& Silk 20031. However, a galactic spectrum contains age- 



sensitive features (e.g. the H/3 absorption line; Worthey & 
'Ot taviani]|1997[) and me tallicity-sensitive features (e.g. the 



MgFe] index; Gonzalez 1993 Thomas et al. 20031, which 



can be useful in breaking this degeneracy. At lower spec- 
tral resolutions the absorption features produced by massive 
stars, for example Balmer lines, can be weakened by coin- 
cident emission within a starburst galaxy. However, higher- 
order Balmer lines are less prone to this effect because the 
strengths of corresponding nebular Balmer emission lines 



strongly decrease with decreasing wavelength ( Gonzalez Del- 
gado ct al. 1998 ). Therefore, a combined analysis of higher- 
order Balmer lines and other tracers may give an indication 
of the ages and metallicities of the dominant underlying stel- 
lar populations in a galaxy. To obtain a range of possible 
ages and metallicities in NGC 4449 we used the spectral 
fitting code STA RLIGHT (v04; |Cid Fernandes et al.|[2005 



Asari et "ar]|2007 1 and the global optical spectrum acquired 



by Kennicutt (19921. This analysis will provide a first indi- 
cation of the star formation history, which will form basis 
for more detailed stellar population fitting in Section [3. 7. 3| 
From the observed optical spectrum, sampled every 2 A, 
we selected the range 3700-4900 A for fitting, with the con- 
taminating Balmer emission lines H/3-H?7, as well as the un- 
resolved [O n]AA3726,9 doublet and [Ne iii]a3869, masked out. 
We adopted 280 kms~^ for the global recession velocity by 
direct measurement from the observed spectrum, but the re- 
cession velocities are generally lower and vary as a function 
of position ( Valdez-Gutierrez et al.|[2002 l. The observed in- 
ternal velocity dispersion, on the other hand, is ~ 30 kms~^ 
( Fuentes-Masip et al. 2000 1 . Since the grating selected by 
Kennicutt ( 1992 1 offered a resolving power of only R ~ 900, 



visited in Section [3.7.1[ Using similar assumptions and the 



the velocity dispersion cannot be resolved by the observa- 
tions and the code cannot take advantage of the more suit- 
able high-resolution stellar spectra for fitting. Nevertheless, 
the distribution of stellar populations resulting from spec- 
tral fitting performed here is not sensitive to the assumed 
internal velocity dispersion. 

We used the synthetic single stellar population (SSP) 
spectra of |Bruzual fc Chariot (2003) with R ~ 2000, which 
were generated with the stellar initial mass function (IMF) 
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Figure 3. One of 12 STARLIGHT fits, in tiie range 3700- 
4900 A, to the global spectrum of NGC 4449 acquired by 
|Kenmcutt| | |1992[ l. The input spectrum (blaclc line) was de- 
composed using 33 different stellar ages at Z/Tiq = 0.2. Con- 
taminating emission lines in the observed spectrum, shown in 
blue, were masked out. See text for more details. 





Mass fraction by age of stellar population 


Metallicity 


Young 
{< 10 Myr) 


Intermediate 
(~ 10^ Myr) 


Old 

(> 10=* Myr) 


0.2Zo 


1% 


20-25% 


60-75% 


O.4Z0 


<2% 


2% 




Zo 


<1% 


<7% 





Table 3. Summary of results from STAR LIGHT fits to th e 
global spectrum of NGC 4449 acquired by |Kennicutt| ( |1992l l. 
The figures show the percentages of the total mass of stars 
contributed by populations in the given age and metallicity 
bin. 



of Chabrier (20031 and provided as part of the STARLIGHT 



package. This assumed IMF is very similar to and consistent 



(1991 



with the IMF of |Kroupa et al.| ( |1991[ |2001[ ) adopted in this 
work. STARLIGHT was run with 12 sets of 5-99 stellar 
templates, or 'bases', to probe the star formation history 
at different sampling resolutions of age and metallicity. The 
lower and upper limits of the stellar ages were 1 Myr and 
12 Gyr and 1-3 metallicities, corresponding to 0.2-1.0 Z©, 
were tested at a time. One of the fits is presented in Fig. |3] 

Our fits showed a significant spread in age-metallicity 
pairs at high resolutions of age and metallicity. Therefore, 
in Table |3] we present a summary of the fitted populations, 
grouped into 'young', 'intermediate' and 'old' and based on 
results from all fits. Table[3]shows that the most metal-poor 
populations are relatively well constrained and represent al- 
most all of the stellar mass in NGC 4449. We find no evi- 
dence for a significant mass of stars with metallicities above 
0.2Zq amongst the old (> 1 Gyr old) population. Based on 
our results, we estimate the first onset of star formation at 
approximately 12 Gyr ago. Because STARLIGHT tended 
to select the oldest available template for the oldest pop- 
ulation, we emphasise that this result is approximate and 
onset epochs ranging from 10 Gyr ago to the age of the 
Universe should be considered equally plausible. Our results 
also consistently point to an age of 400 ± 100 Myr as being 
representative for the intermediate population. 

These results are supported by the resolved stellar pop- 



ulation study of Annibali et al. ( 2008 1, suggesting continuous 
star formation activity for at least 1 Gyr, and by the results 
of Bothun (19861, who estimated that the underlying old 



stellar component has a mean age of 3-5 Gyr. A very young 
population of stars, on the other hand, is inferred from the 
presence of Wolf-Rayet stars ( Martin fc Kennicutt|1997 1 . 

Therefore we assume that star formation started not 
later than ~4 Gyr ago and continued until very recently 
or continues to the present day. We also assume three con- 
tinuous star formation episodes. This assumption reduces 
the complexity of the model, and, at the same time, hints 
at a possible star formation history. The continuous-episode 
approach can be viewed as a generalisation of a starburst 
approach, which is discussed in Section [4. 1| 



3.5 Interstellar dust 

3.5.1 Distribution and composition 

The dust distribution is assumed to follow that of the gas, 
with the dust to gas mass ratio (abbreviated here as DGR or 
Afdust/Afgas) initially kept constant throughout the galaxy 
(however see, e.g.|Bianchi|2008, Baes et al.|2010 MacLach 



Ian et al.|2011||Schechtman-Rook et al.|2012[|de Looze et al 



2012 1. We assumed two dominant dust species: amorphous 
carbon ( |Hanner||1988[ ) and silicates ( |Laor fc Dra inc 1993|. 
To limit the number of free variables in the model, we further 
assume a constant mass ratio of carbonaceous dust to sili- 
caceous dust (hereafter referred to as the carbon-to-silicate 
ratio) as a function of position. In what follows, we intend 
to illustrate explicitly the evolution of the carbon-to-silicate 
ratio for a galaxy continuously forming stars to allow us to 
fix this ratio for NGC 4449. We refer the reader to IDwekl 



( |1998|), [Morgan fc Edmunds] (|2003| ), [Galliano et al.| ( |2008a l 
and [Dwek fc Cherchneff[ ( [2011[ ) for more complete studies of 
dust evolution in the ISM. 

The refractory dust in the ISM is believed to form 
in the outflows from asymptotic giant branch stars (AGB 
stars; e.g. Matsuura et al. 2009 for the LMC), in type II 



core-collapse supernova (CCSN) ejecta (e.g. a multi-epoch 
study of SN 2004et by [Fabbri et al.|2011|), and has al s o been 
suggested to form in molecular clouds ( Draine|2009 1. Dwek 



( 1998 1 has suggested that low-mass AGB stars are the main 
source of carbon dust, while type II SNe are the main source 
of silicate dust. Thus, the relative abundance of carbon and 
silicate dust depends on the dynamics of dust formation by 
AGB stars and supernovae, and on dust destruction pro- 
cesses and thus on the lifetimes of individual dust species. 

Our calculations assume the dust yields of [Ventura et al.| 
( |2012[ ) for AGB stars, a constant and continuous star forma- 
tion activity (Section [3.4| above), a constant metallicity, the 
initial mass function (IMF) of [Kroupa et al.[ ( [1991[[2001[ ) for 
stellar masses 1.5 M© < M < 40 Mq, and ignore dust de- 
struction processes including consumption by ongoing star 
formation. 

It is informative to first consider a flat distribution of 
yields across all SN masses and the carbon-to-silicate ratio 
inferred for SN 1987A by [Matsuura et al.[ ((20TT1 Model 



1). Fig. |4] shows the carbon-to-silicate ratio for a galaxy ac- 
tively forming stars for a continuous period of 4 Gyr. Since 
the data of [Ventura et al. (20121 do not include AGB stars 
of mass lower than 1.5 Mq, our models cannot be reliably 
extrapolated beyond the first 4 Gyr. At early times, the 
carbon-to-silicate ratio is fixed by the assumed dust compo- 
sition of SN ejecta. When the SN dust dominates the global 



© 2002 RAS, MNRAS 000,[lH2T] 



8 O. L. Karczewski et al. 



dust budget, the ratio is not expected to change significantly 
with time (Fig. |4] left, red solid line). If AGB dust domi- 
nates dust production, the carbon-to-silicate ratio initially 
decreases after ~ 40 Myr due to the predicted injection of 
silicates by super- AGB stars. After ~ 400 Myr, less mas- 
sive AGB stars begin injecting carbonaceous dust leading to 
carbon-dominated ISM dust after ~ 1 Gyr (Fig. [4] left, blue 
dotted line). 

A more realistic model includes variations in dust mass 
production as a function of CCSN progenitor mass, based on 
CCSN element yields o f|Woosley fc Weaver] ( |1995[ ). Since the 



metallicity assumed by Ventura et al. ( 2012 1 is Z/7iq = 0.05, 



we chose the closest Z/7iq = 0.1 Ni-producing models 
from [Woosley fc Weaver] ( |1995 !) for all stellar progenitor 



masses, thus ignoring fallback of heavy elements (Moriya 
et al.||2010| ). We further scaled the CCSN elemental yields 



by a constant factor of 0.2 to bring the predicted dust masses 
into approximate agreement with the dust masses inferred 



from a range of supernovae (|Wesson et al. 2010 


Temim 


et al.|2006| 


Barlow et al.|2010||Matsuura et al..2011| 


Gomez 


et al. 


2012 


1. This scaling factor may be physically inter- 



preted as the adopted condensation efficiency of dust. As 
amorphous olivine is observed to dominate the silicate mass 
in the diffuse ISM ( Kemper et al.|2004 |, we assumed silicate 
stoichiometry of MgFeSi04 for calculations of silicate dust 
masses. 

The results in Fig. [4] (right) show that the carbon-to- 
silicate ratio increases with time, owing to more carbon-rich 
dust being produced by lower mass SNe. This is strength- 
ened at later times by carbon-rich dust contribution from 
the less massive AGB stars. 

In this model a single metallicity of Z/Z0~O.l was 
assumed for the dust. However, the yields of [Woosley fc| 
Weaver ( 1995 1 suggest that as metallicity increases, the rel- 



ative production of carbon by CCSNe also increases. The 
opposite trend is suggested for AGB stars. In low metal- 
licity environments AGB stars are expected to produce less 
silicates, while maintaining the same or an enhanced produc- 
tion of carbon dust ( Sloan et al.|2006 20121 than in higher 
metallicity environments. Since these trends operate on dif- 
ferent time-scales, their combined effect is likely to result in 
a carbon-to-silicate profile which is more steeply increasing 
than the one shown in Fig. [4] 

Therefore, given our assumptions about metallicity, the 
limited period of modelled star formation activity, and the 
uncertainty of the degree of dust enrichment by SNe, our 
final carbon-to-silicate mass ratio of 0.23 (Fig. |4] right) is 
likely to represent the lower limit of relative carbon content 
in a continuous star formation scenario. NGC 4449, with 
log(0/II) -I- 12 = 8.23, is similar in terms of metallicity to 
the LMC ( [Russell fc Dopita|1992[ ), whose dust content sug- 
gests a slightly higher carbon-to-silicate mass ratio of 0.33 



(Weingartner & Draine 20011. If our model is allowed to 
evolve beyond the first 4 Gyr, the relative carbon content in 
NGC 4449 would be expected to match or exceed the LMC 
value, depending on the age of NGC 4449. This observa- 
tion may warrant the use of metallicity-based assumptions 
about the carbon-to-silicate ratio in future studies. There- 
fore, given that our spectral fitting results discussed in Sec- 
tion [3]4] point to an early onset of star formation, we adopt 
the LMC carbon-to-silicate mass ratio of 0.33. 



3.5.2 Grain size distribution 

While the absence of distinct silicate emission features in 



the Spitzer Infrared Spectrograph (IRS) spectra (Paper III 
may indicate a small abundance of Si, our chemical evolution 



models (Section 3.5.1 above) predict a non-negligible mass 
of silicate dust. If all dust is assumed to reside in a thin layer 
around the central ionizing source, this would impose an up- 
per limit of 20 per cent on t he mass fraction of Si in 'ultra- 
small' silicate grains (< 15 A; |Li fc Draine|200l| ) and suggest 



that the contribution of hot silicate dust to the global dust 
budget is relatively small (e.g. Smith et al.|2010 l. However, 
NGC 4449 is a face-on irregular galaxy and it is reasonable 
to assume, to a first-order approximation, that the radial 
distribution of dust follows the radial distribution of gas. 
Therefore, the absence of distinct silicate emission features 
in the global IRS spectrum may indicate effects resulting 
from geometry, for example extinguishing the emission from 
small grains by self-shielding. 

We assum e that the dust grain size distribution follows 



n(a) cx a (Mathis et al. 19771 and we adopt the same 



range of grain sizes a between 0.005 /j,m and 0.25 /xm for 
both amorphous carbon and silicates. 



3.5.3 Poly cyclic aromatic hydrocarbons 

Prominent emission features in the mid-infrared, also known 
as the unidentified infrared bands (UIBs), are visible in the 
spectra of many types of objects, and are also detected in 
the Spitzer spectrum of NGC 4449 ( PaperTT] ). They are 
believed to originate from the C-H and C-C vibrations of 
large organic molecules, or polycyclic aromatic hydrocar- 



bons (PAHs; e.g. AUamandola et al. 1985 Peeters et al 



2002 Bauschlicher et al. 20091. The relative strengths of 



these emission features are believed to refiect the particle 
size distribution and the local physical environment of the 



PAH molecules (e.g. Galliano et al. 2008b Bauschlicher et al. 
2009] ). 

We assumed the mathematical description of pure ion- 



ized PAHs from Draine fc Li ( 2007 1 with the grain size dis- 



tribution of Weingartner & Draine 



(20011 for grain sizes of 



3.5-30 A. The PAH masses in this approach are not con- 
strained by abundances (cf. Zubko et al.||200"4 i. Since MO- 
CASSIN currently allows only one grain size distribution per 
simulation, the PAHs have to be modelled separately. 



3.6 Numerical setup 

A series of preliminary simulations showed that the best 
compromise between high spatial resolution and reasonable 
computing time is obtained for a grid of 80 x 80 x 80 cells, cor- 
responding to 0.5 million grid cells. In a spherically symmet- 
ric setup, MOCASSIN uses a Cartesian x, y, z grid, with cells 
populating one octant of a sphere resulting in an eight-fold 
reduction of computing time for a given spatial resolution. 
The grid cells were not equal in size, and the resolutions used 
were 15 pc, 50 pc and 185 pc for the inner 0.5 kpc region, 
intermediate radii up to 2.5 kpc, and for outer radii, respec- 
tively. The number of grid cells and the resolutions ensure 



adequate sampling of the gas density profile (Section 3.2 I as 
a function of radius. A single ionizing source, describing the 
stellar content of the entire galaxy, was placed at the origin 
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Figure 4. A simple model of the evolution of the global relative carbon to silicate dust content of NGC 4449, assuming continuous 
star formation and a constant metallicity Z/Zq ~ 0.1. The left panel shows the evolution of the carbon-to-silicate ratio for 
supernova dust yields identical to those derived by |Matsuura et al.| | |20lT| | for SN 1987A (solid line), adopting their lower limit 
of -Mjjust = 0-4 Mq and a carbon-to-silicate ratio of 0.21 (their Model 1). Similar models with dust masses scaled to 0.01 Mq 
(dashed line) and 0.001 Mq (dotted line) are also presented, where the 0.001 Mq model may be viewed as one where core-collapse 
supernovae are not important in global dust production. The right panel shows the evolution of the carbon-to-silicate ratio for 
supernova dust yields based on the elemental yields of |Woosley fc Weaver, ( |1995| solid line). The thin solid and dashed lines show, 
respectively, the relative masses of carbon and silicate dust produced. 



of the coordinate system. The physical inner and outer radii 
of the galaxy were defined as 1 pc and 3300 pc. 

The computation time for a single gas and dust simu- 
lation resulting from the setup described in this paper was 
~ 9 h using 8 CPUs with 8 GB of RAM per CPU. 



3.7 Variables and convergence criteria 

As described in Section |3.1| and shown in Fig. [2j the three 
variables, the global bolometric luminosity L*, the DGR and 
the assumed star formation history (SFH), are iteratively 
adjusted until all obscrvables are well matched. Of these, 
Li, is the most independent and can be fixed at early stages 
of the modelling. 

The DGR and the SFH are not independent, as a higher 
dust content increases the extinction at the UV and optical 
wavelengths and therefore requires more ionizing photons to 
maintain the observed continuum levels. The SFH in itself 
is a multi-dimensional function, which can be explored to 
reach the most probable solution. We use the general pic- 
ture of the SFH emerging from Section |3.4| as a constraint 
for constructing trial three-episode star formation scenarios. 
The steps described below may be viewed as an attempt to 
minimise a multi-dimensional mathematical expression by 
examining its partial derivatives locally. 

3. 7. 1 The bolometric stellar luminosity L* 

A large fraction of all the energy emitted by a starburst 
galaxy comes from the youngest and most massive stars. 
These stars contribute to the observed UV and optical con- 
tinuum and govern the observed global ionization states of 
species present in the ISM. The same radiation field is also 
partially absorbed by the interstellar dust and re-radiated 
thermally in the IR. The degree of ionization of nebular 
species and the degree of reprocessing of radiation are also 
dependent on the dumpiness of the ISM. 

Therefore, by focusing initially on the youngest stellar 



component and considering a representative set of star for- 
mation scenarios it is possible not only (i)to further simplify 
the SFH by constraining this youngest component, but also 
(ii) to fix the filling factor e, (Hi) to fix L* and (iv) to pro- 
vide an initial estimate for the DGR. Two main constraints, 
namely, the absolute level of the UV continuum and the 're- 
cent' star formation rate (SFR) as traced by Ha or a similar 
tracer, must both be matched by the total stellar luminos- 
ity and the average age of the youngest population given the 
adopted filling factor. 

The 'recent' star form ation rate for NGC 4449 is es- 
timated at ~0.5 Moyr"^ |Hunter et al.||l986l |l999j) for a 



Salpeter IMF ( Salpeter 1955 see discussion in Section 3.7.2 1 
Since [Popita et al ( 2006| showed that for a SSP all ionizing 
photons are emitted within 10 Myr, we adopt 10 Myr, the 
lifetime of O- and B-type stars with masses > 13 Mq, as a 
representative period of continuous star formation for the 
youngest stars. 

The trial star formation histories (here also referred to 
as 'scenarios') were generated using STARBURST 99 (|Lei- 



therer et al. 19991 v6.0.2) with the IMF of Kroupa ( |Kroupa 
et al.||199ir200l|, along with t he Pauldrach/Hillier model 
atmosp heres ([Smith et al.|200 2), the updated Padova AGB 



tracks ( Vazquez fc Leitherer|,2005) and an assumed single 
constant metallicity of Z/Zq = 0.4 (cf. Table [3|. 

For a representative set of scenarios consisting of three 
star formation episodes, we first varied the degree of dumpi- 
ness, represented by e, to allow enough gas to be ionized and 



produce the observed nebular emission line intensities ( Kob- 



ulnicky et al. 1999 1 to within a factor of a few. We found 



that e = 0.033 is most representative, which corresponds 
to an average nc of 41 cm""^ within the modelled ionized 
region. This choice is discussed further in Section [Xs] 

We constructed a grid in parameter space formed by a 
representative set of scenarios consisting of three star forma- 
tion episodes, a range of SFRs of the youngest population 
which is continuously forming stars for a period of 10 Myr, 
and a small range of initial values for and the DGR. 
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The predictions were then compared with observations. In 
particular, 

1. the predicted UV continuum level was compared with 
the observed SED (Table [2|, 

2. the estimated total number of ionizing photons, QiYiP), 
was compared with that derived from the integrated Ha 
luminosity I Hunter et al.|1999 l and 

3. the predicted nebular emission line fluxes were com- 



agreement with the distance-cor rected and IMF-cor rected 
estimates of 0.22-0.33 M pyr"^ ( [Hunter et aL||l986| , 0.30 
M0yr"^ ([Hunter et aI]|T999), and 0.21 Mpyr^^ (based on 



pared with the measurements of Kobulnicky et al. ( 1999 1 



The best-fitting models were found by minimising resid- 
uals between the model and the observations in a way similar 
to minimisation. At this stage of modelling more weight 
was given to line intensities and to the level and shape of the 
continuum in the UV. To reduce the extent of the parameter 
space and the overall computation time, the iterative scheme 
in Fig.|2|was not allowed to advance automatically. Instead, 
after each iteration, candidate parameter sets were inspected 
to identify trends in each variable. The most promising pa- 
rameter combinations were then manually expanded into a 
higher-resolution parameter space and iterated according to 
Fig. [2] 

Our results suggest an on-going star formation activ- 
ity with recent SFR 0.28 M0yr"\ 5.7 x 10^ Lq 
and an initial DGR ~ 1/250, thus allowing us to fix the 
youngest population and Li,. The empirical L*, computed 
by integrating measurements in Table [2| between 0.15 /im 
and 3.5 /xm, is 3.1 x 10^ Lq. This suggests that a significant 
fraction of the total luminosity originates from the recently 
formed massive stars. 



3.7.2 The initial mass function 

The recent SFR of 0.28 MQyr~^, suggested by our initial 
models, is in good agreement with previous estimates (see 



below), when the revised distance of 3.8 Mpc ( Annibali et al 
[2008[ ) and differences in the assumed IMFs are taken into 
account. 

The Salpeter IMF, often used in the formulations of 
the SFR (e.g. |Kennicutt"1998|), is defined as iV(M) dM oc 



;.IKennicutt"1998|), 
((Saipctcr 1955|7^I 



M'^--''^ AM ([Saipctcr 19551). However, the Kroupa IMF, 
which is assumed in this work, is defined in two intervals, 
with a flatter distribution of the number of stars at lower 
stellar masses. Therefore, the Salpeter IMF will overestimate 
the number of lower-mass stars compared to the Kroupa 
IMF. This will in turn overestimate the total mass of stars 
resulting in a higher prediction of the SFR. 

We generated distributions of stellar masses normalised 
to give the same number of stars at a stellar mass of 25 
Mq, which we assumed to be representative of the popula- 
tion generating the ionizing photons and giving rise to the 
observed Ha emission. We computed the total mass of stars 
from N{M)M AM for 0.1 ^ M/Mq ^ 100, which yielded 
a factor of 1.54 difference between the two IMFs. Assum- 
ing that the integrated luminosity of a galaxy scales linearly 
with the total mass, the SFRs based on the Salpeter IMF 
are therefore likely to be a factor of ~ 1.5 higher than if the 
Kroupa IMF were used. This result is also evident from the 
results of Dwek et al. (20111, who tabulate the masses of 



high ionization potential neon lines; Paper II I 



all stars born per one supernova event for a range of stellar 
IMFs. 

Overall, our recent SFR of 0.28 M0yr~^ is in good 



3.7.3 Three episodes of star formation 

To constrain and the youngest stellar population, we as- 
sumed a representative older stellar component. Since the 
contribution of stars older than ~ 10 Myr to the ionization 
structure, to the number of ionizing photons, or to the far- 
UV continuum, is very small, the older populations can be 
studied in more detail without affecting the validity of our 
findings from Section [3.7.1| 

The shape of the observed SED from the UV to the 
near-infrared (NIR) can be fitted by varying the lengths and 
characteristic epochs determining the remaining episodes of 
continuous star formation. Initially, we tested a set of 36 
scenarios, which consisted of all combinations of: 

1. three epochs corresponding to the onset of star forma- 
tion (4, 8 and 12 Gyr ago), 

2. four epochs corresponding to the transition between 
the first (old) episode and the second (intermediate) episode 
of star formation (100, 200, 400 and 1000 Myr ago), 

3. three star formation rates ranging from 0.05 M0 yr~^ 
to 0.15 M0 yr^^ further defining the second episode. 

The star formation activity during the third (young) episode 
remained fixed in both duration and magnitude, and L-,, 
represented the luminosity of all stellar components, as de- 
scribed in Section [3.7.1[ Therefore, the star formation rate 
of the first (oldest) episode was not an independent variable, 
but depended on and the durations and star formation 
rates in items 1-3, for each scenario. 

The resulting low-resolution SEDs were compared with 
the observed SED to select the best matches. Afterwards, 
following the scheme shown in Fig. [2] the transition epochs 
and the star formation rates inferred from the best matches 
were expanded into a narrower, but higher-resolution pa- 
rameter space, analogous to conditions in items 2-3. 

The criteria for finding the best-fitting combinations of 
stellar populations were similar to those used to find L*, 
where we minimised residuals between the model and the 
observations. At this stage, more weight was given to the 
number of ionizing photons (3(H°), the line intensities and 
the level and shape of the entire stellar continuum from the 
UV to the NIR. As before, to reduce the extent of the pa- 
rameter space and the overall computation time, the itera- 
tive scheme in Fig. [2| was not allowed to advance automati- 
cally. Instead, after each iteration, candidate parameter sets 
were inspected (i) to identify trends in each variable and 
(a) to verify if the fit could be improved while still satis- 
fying the observed number of ionizing photons Q(H''). The 
most promising parameter combinations were then manually 
expanded into a higher-resolution parameter space. 



3.7.4 ^ two-zone solution for Mdust/M^^s and modelling 
PAHs 

NGC 4449 is modelled as a centrally-concentrated, 
spherically-symmetric galaxy, even though its morphology 
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is highly irregular (cf. Fig.[T]). Therefore, spherical symme- 
try inevitably misrepresents the interstellar radiation field 
(ISRF) within the galaxy. The single central source is much 
more luminous than any individual source in the real galaxy, 
as it must deliver enough energy to the correct mass of gas 
and dust to reproduce the observed strengths of the emission 
lines and the observed SED. Consequently, the conditions in 
the immediate proximity of the single central source result 
in unrealistically high ionization of the gaseous species and 
vaporisation of any existing dust grains. This effect is dis- 
cussed in more detail in Section [4.31 

In order to fit the dust emission in the FIR, we divided 
the galaxy into two zones: an inner zone with a lower DGR 
and an outer zone with a higher DGR, which ensures realis- 
tic absolute masses of hot dust near the centre. For techni- 
cal reasons this was preferred to an alternative approach, in 
which the radial distribution of dust density is defined di- 
rectly through the MOCASSIN input parameter Ndust. We 
found that for larger systems, such as galaxies, defining dust 
distribution through MdMg gives best performance. 

The best-fitting models were repeated with pure ion- 
ized PAH dust and the PAH grain size distribution (Sec- 
3.5.3|l, and compared with the observed global SED 



tion 



(Table [2|) and the global Spitzer /IWSt spectrum (Paper III 
to estimate the total mass of PAH molecules. In the final 
step, the resulting spectra were combined with the corre- 
sponding amorphous carbon and silicate models. 



4 RESULTS AND DISCUSSION 
4.1 Preliminary models 

Our preliminary models included both starburst-only sce- 
narios and continuous star formation scenarios. In the 
starburst-only scenario, a fit satisfying all observational con- 
straints can be obtained for three representative starbursts 
occurring 3 Myr ago, 100 Myr ago and 4 Gyr ago, con- 
tributing to the total stellar mass in a ratio 1:200:1000. 
These models were later generalised by replacing starbursts 
with longer star formation episodes. For a scenario assuming 
two continuous star formation episodes the inferred star for- 
mation rates were ~0.25 Mq yr~^ and ~0.09 M© yr~^ be- 
tween ~ 6 Gyr ago and ~ 120 Myr ago, and ~ 120 Myr ago 
and the present day, respectively. The star formation his- 
tory in these models was constrained only by the observed 
SED and by the observed emission line intensities. We note 
that the representative populations and star formation rates 
are not too dissimilar to the assumptions described in Sec- 
and the recent star formation rate of 0.28 M© yr~^ 



3.4 



tion^ 

(Section 3.7.2[ |. Therefore, it may be argued that combin- 
ing the three convergence criteria, namely, (i) matching the 
observed SED, (ii) matching the total observed rate of ioniz- 
ing photons, Q(H*'), and (Hi) matching the observed nebular 
emission line fluxes, may give informative results, even for a 
relatively simplified model. 



4.2 Final model 

Fig. [5] shows the best-fitting MOCASSIN model of 
NGC 4449 satisfying the observational constraints over the 



entire wavelength range from the UV to sub-mm. The in- 
dividual SEDs arising from the assumed three continuous 
star formation episodes are also shown for comparison. The 
input parameters and the results are summarised in Table [4] 
and the details of the best-fitting parameters describing the 
star formation history are given in Table |5] 

The three episodes of star formation shown in Fig. [5] 
suggest that most of the UV emission, and most of the 
processed radiation emitted in the FIR, originates from the 
intermediate-age stars of ages between 10 and ~ 400 Myr. 
The stars produced in the first episode, i.e., those older than 
~ 400 Myr, dominate in the optical and in the NIR. The 
youngest stars of ages less than 10 Myr contribute mainly 
in the far-UV. 

Additionally, Fig. [6] shows the effects of varying the in- 
put parameters for the same best-fitting model (red solid 
line). For example, the effect of varying L* is a vertical shift 
of the predicted SED in the stellar part (Fig. [6] top panel) 
and a corresponding change in Q(H''), as the IMF- weighted 
distribution of stellar luminosities is scaled by a constant. 

The predicted and observed global emission line intensi- 
ties, relative to H/3, are shown in Table[6] The discrepancies 
between the predicted and the observed nebular line inten- 
sities of some species, most notably [S ll] , result from the 
simplifying assumptions about the geometry and the degree 
of dumpiness, and are discussed in Section [4. 3[ Because the 
current version of MOCASSIN offers no treatment of PDRs, 
the PDR line [O l] in Table |6] is underpredicted. However, 
MOCASSIN is in the process of being expanded to accom- 
modate PDRs via the recently developed three-dimensional 
PDR code 3D-PDR ( [Bisbas et a"L||2012 K 



4-2.1 Stellar populations 

In our three-episode model a stellar mass of ^ 1 X 10^ Mo 
is produced at rates 0.25-0.10, 0.14 and 0.28 M© yr~^ over 
three periods, as shown in Table |5] The best-fitting SED in 
Fig. [5] shows that this model does not fully account for the 
emission between 7000 A and 2 /im, which may be a result of 
the assumption of only three episodes of star formation with 
a constant SFR for the duration of each episode. However, 
we note that a non-trivial three-dimensional dust distribu- 
tion, where the degree of obscuration does not follow the 
distribution of gas or stars, may have a significant effect 
on the observed attenuation (e.g. Witt et al.||1992 Baes & 
|Dejonghe|2"000l ). 

We found that the results are degenerate for the oldest 
stars and therefore the first onset of star formation 4, 8 
or 12 Gyr ago is equally plausible, resulting in a range of 
predicted star formation rates. The spectral fits performed 
with STARLIGHT and discussed in Section [3.4| suggest an 
onset 12 Gyr ago or earlier. For an onset of star formation 
12 Gyr ago, the average star formation rate for the duration 
of the first episode is 0.09 M0yr~^. Because the mass of 
recently formed stars is small in comparison (cf. Tables [3] 
and|5|, it is possible that up to 95 per cent of the total stellar 
mass Mi, in NGC 4449 was produced at higher instantaneous 
rates in an initial starburst. 

The predicted onset of the second episode, 300-400 Myr 
ago, coincides with a possible encounter with the galaxy 
DDO 125 400-600 Myr ago | |Theis fc Kohle|2001| ). Given the 
simplifications of our model it would be difficult to quantify 
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Parameter 



Value 



Comments 



geometry 
grid size 

number of photons 



physical radius 
gas distribution 

Mgas 

filling factor (e) 
chemical elements 
Q(H«) 



dust composition 
dust grain sizes 
PAH grain sizes 



Li, 

MpAH 
A^dust/Mgas 



Numerical setup 

3D, spherically-symmetric 
80 X 80 X 80 
2 X 10® 



ionizing source at (x, y, z) = (0, 0, 0) 
over all wavelengths 



Observationally constrained parameters 

3.3 kpc 

exponential-like 



0.55 X 10^ M© 
0.033 

H, He, C, N, O, Ne, S 
4.84 X 10^2 s-i 



derived from the H I profile of 
Swaters et al.| l |2002[ | 

computed for r 3.3 kpc; molecular 
gas not included 

abundances constant throughout galaxy 
based on Lhq l |Hunter et al.|1999| 



Adopted dust characteristics 

amorphous carbon and silicates (1:3) 
0.005-0.25 fim (20 sizes) 
3.5-30 A (10 sizes) 



Hanner| ||1988[l;|L aor fc Draine| l |l993 I 

athis et al.'(197^ 
FAHs modelled separately; 



Weingartner & Draine ( 2001 1 



Best-fitting parameters 

5.7 X 10^ L0 

1.05 ±0.15 X 10^ Mq 

2.9 ±0.5 X lO'' Mq 

0.058 ± 0.005 X 10^ Mq 

1/190 



Draine &; Li 



including PAH grains 
ionized PAH grains; 
effective; 
1/1000 in centre, 1/182 otherwise 



1 2007 1 



Table 4. Summary of the input parameters and the results from the best-fitting MOCASSIN model of NGC 4449. 
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Figure 5. SED of the best-fitting MOCASSIN model of NGC 4449. The three star formation episodes were modelled simultaneously, 
but their individual SEDs are also shown for comparison. The photometric measurements have been corrected for foreground 
extinction, and the mid-IR spectrum acquired with Spitzer /TRS is discussed in detail in [Paper II| 
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the effects of this postulated encounter on star formation 
rates. However, we note that our STARLIGHT fits (Sec- 
tion 3.4 1 select a population ~ 400 Myr old as being rep- 
resentative for the intermediate-age populations. Similarly, 
the three-episode MOCASSIN fits to the observed stellar 
SED converge to 300-400 Myr ago as the transition epoch 
between the first episode and the second. 

Our model suggests that in the last 10 Myr star for- 
mation has been taking place continuously at an average 
rate of 0.28 M q yr"^ (0.42 M© yr"^ assuming Salpeter IMF; 
Section |3.7.2[ ). This recent SFR, obtained by fitting both 
the UV continuum and the emission lines (Section 3.7.1 1, is 



consistent with other estimates based on the UV continuum 
emission (0.44 M© yr~^; using data in Table [2] and following 
the formalism of Kennicutt|1998j) or em ission line intensities 
i Hunter et ar||1986[ [Hunter et al.|1999[ [Paper \\\ . 

The FIR emission arises from dust heated by stars and 
can also be used to estimate the recent SFR. Integrating 
the best-fitting SED in Fig. [5| between 8 and 1000 /xm yields 
/(FIR) 3.6 X 10"^^ Wm"^. The corresponding recent 
star formation rate is 0.28 M0yr~^ assuming a Salpeter 
IMF ( Kennicutt|1998 1, which is significantly lower than the 
estimate of 0.44 M0yr~^ based on the UV emission. In- 
deed, the ratio SFR(FIR)/SFR(UV) for NGC 4449 based 
on the data in Table [2[ and on the fit presented in Fig.[5j is 
0.28/0.44 = 0.64, suggesting that approximately 35 per cent 
of the UV radiation does not contribute to the heating of 
dust. This UV 'leakage' (e.g. Relano et al.|2012 l is observed 
also in Haro 11 and NGC 4214 ( [Cormier et al.[r2012[ [Her- 
|melo et al.[[2013| ), and may result from porosity of the ISM 
or different distributions for the stars and dust. 

The low ratio SFR(FIR)/SFR(UV) may also suggest 
that, globally, in NGC 4449 the dust is heated by younger 
stellar populations. Similarly, a spatially-resolved study by 



Galametz et al. (20101 suggests that the distribution of 



cooler dust within a different dwarf galaxy, NGC 6822, may 
be correlated with star formation activity. In spiral galaxies, 
on the other hand, several studies suggest that the dust is 
heated by both the younger and the evolved stellar popu- 
lations (e.g. Bendo et al.[[20l"0 Boquien et al.[[2011 Bendo 



et al.|2012a] 



4-2.2 Interstellar dust 

We infer a dust mass of 2.9±0.5 x 10® Mq by making assump- 
tions about the distribution of dust and its characteristics: 
the carbon-to-silicate ratio and the grain size distribution 
(of. Section [3^ and Table jif. 

As explained in Section |3.7.4[ in modelling the dust 
emission within NGC 4449 we adopted two zones with differ- 
ing DGRs. Fig. [6] (middle panel) shows the two components, 
modelled simultaneously, in green and blue. The combined 
effect of the best-fitting DGRs and the adopted distribution 
of gas (Section 3.2 \ is a relatively constant distribution of 
dust as a function of radius, varying by less than a factor 
of three between the centre and r = 3.3 kpc, and yield- 
ing Mgas/190 within the modelled region. The need for two 
zones demonstrates that the actual distribution of dust is 
significantly different from the distribution of gas, even for 
one-dimensional azimuthally-averaged profiles. Fig. [6| also 
shows results for a smaller and a larger overall DGR, sug- 



gesting that the uncertainty in the best-fitting dust mass is 
~0.5 X 10** Mq. 

The composition of dust has a major infiuence on the 
amount of attenuation produced at UV and optical wave- 
lengths, and therefore on the total amount of dust required 
to reproduce the observed FIR emission. Fig. [5| shows that 
the photometric measurements agree with the model near 
the prominent absorption feature at 2175 A, which is at- 



tributed to carbon dust (e.g. Fitzpatrick & Massa 1986 1 



The Spitzer /IRS spectrum, on the other hand, is domi- 
nated by the PAH emission features and shows no discernible 



10 /xm silicate feature ( Paper II I, and no additional underly- 



ing component at 10 /xm was needed to fit the PAH emission. 
Therefore, while appreciating that the assumed carbon-to- 
silicate ratio of 1:3 may be underestimated (also see discus- 



sion in Section 3.5 I, we conclude that the adopted ratio is 



consistent with observations. We find a contribution from 
PAHs to the total mass of dust of 2 per cent, which is con- 
sistent with the MpAn/Afdust ratios expected for a galaxy 
with 1/3Zq ( [Galliano et al.|2008a[ ). 

For comparison, in Fig. [6[ (bottom panel, dotted line) 
we show the SEDs corresponding to different compositions of 
dust for the same total mass of dust. Additional carbon-rich 
dust (dotted line and dashed line) produces more emission 
in the FIR. If such carbon-rich compositions were used in 
fitting the SED, the overall required mass of dust would be 
lower. More silicate-rich dust (dot-dashed line), on the other 
hand, would increase the required mass of dust. Dust mass 



estimates in the literature range from 2.0 x 10 Mq Engel 



bracht et al. 20081 to 3.8 x 10" M© (Bottner et al 



20031, 



refiecting different properties of dust adopted in interpreting 
the FIR emission. 

The metallicity of NGC 4449 is similar to that of 



NGC 1705 and NGC 6822 studied previously (O'Halloran 



et al. 2010 Galametz et al. 20101. However, the DGRs, in 



the range 1/80-1/186, derived for these galaxies are higher. 
The higher ratios may result from the assumption of pure 
graphite or pure amorphous carbon dust, and from the 
adopted gas and dust budget. 

We note that the DGR of 1/190 was computed using 
the total derived dust mass for NGC 4449, while the model 
covered the inner radius of 3.3 kpc, the extent of the galaxy 
in the FIR, and as such took into account only about 25 
per cent of the total observed gas (Section [3.2[ ). If we con- 
sidered only the gas present within the modelled region, that 
would require all metals in NGC 4449 to be locked into dust. 
However, NGC 4449 is a highly dynamic system, postulated 



to have been involved in a merger (Theis & Kohle 20011 



and possessing a recently discovered companion ( Rich et al 
[2012[ ). Since its optical body is not likely to be an isolated 
system, it is informative to include the extended H I envelope 
in the calculations, which yields a lower DGR of ~ 1/760. 
We also note that our gas masses do not include the mass 
of H2 due to large uncertainties in the CO-to-H2 conver- 
sion factors for low metallicity galaxies. Inclusion of an H2 
mass of ~ 1.7 x 10* Mq could decrease the global DGR even 
further to ~ 1/820. 

On the other hand, an elevated flux at 850 pm, visi- 
ble in Fig. [5] and not matched by the model, may suggest 
the presence of a significant mass of cold dust (e.g. 'Galliano] 
et al.||2003) . Correctly accounting for the sub- mm emission 
is likely to increase the total mass of dust, and consequently 
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Episode 1 


Episode 2 


Episode 3 




(Old) 


(Intermediate) 


(Young) 


onset / yr ago 


4-12 X 10'' 


300-400 X 10^ 


10 X 10^ 


M* / Mq 


~109 


~ 50 X 10^ 


~3 X lO'' 


SFR / M0 yr-i 


0.25-0.10 


0.14 


0.28 



Table 5. The best-fitting three-episode star formation history 
of NGC 4449 assuming Kroupa IMF. The assumed star forma- 
tion activity is continuous, and the onsets of episode 2 and 3 
coincide with the end of the preceding episodes. 



Line 


Predicted 


Observed 


Ratio 


References 




[%/3 = 1] 


[Inp = 1] 






[O II] A3727 


3.198 


3.891 


0.822 


(1) 


[Ne III] A3868 


0.173 


0.189 


0.916 


(1) 


H<5 A4101 


0.260 


0.279 


0.932 


(1) 


H7 A4340 


0.469 


0.492 


0.953 


(1) 


A4861 


1.000 


1.000 


1.000 


(1) 


[O III] A4959 


0.699 


0.689 


1.015 


(1) 


[O III] A5007 


2.086 


2.069 


1.008 


(1) 


He I A5876 


0.079 


0.079 


1.000 


(1) 


[S III] A6312 


0.008 


0.019 


0.423 


(1) 


[N 11] A6548 


0.102 


0.115 


0.892 


(1) 


Ha A6563 


2.922 


2.866 


1.020 


(1) 


[N 11] A6584 


0.313 


0.338 


0.927 


(1) 


He I A6678 


0.022 


0.028 


0.798 


(1) 


[S 11] A6716 


1.117 


0.476 


2.347 


(1) 


[s 11] A6731 


0.784 


0.334 


2.349 


(1) 


[s iv] 10.5 Atm 


0.264 


0.219 


1.205 


(2) 


[Ne III] 15.6 Aim 


0.228 


0.387 


0.588 


(2) 


[O i] 63 fim 


0.012 


0.413 


0.030 


(2) 


[O III] 88 ^lm 


0.654 


0.647 


1.011 


(2) 


[C 11] 158 ^lm 


1.112 


1.307 


0.851 


(2) 



Table 6. Predicted and observed line strengths relative to H/3 
and the predicted-to-observed ratio for the best-fitting MO- 
CASSIN m odel of NG C 4449. References: (l) |Kobulnicky et al.| 
([T999J; (2) [Paper Il| 



the DGR, and can be crucial for accurate dust mass deter- 
minations (e.g. Galametz et al.|2011 |. 

Finally, the dust mass predicted by our model can be 
used to infer the fraction of metals locked into dust in 
NGC 4449. Adopting the total Hi mass of 2.2 x lO" Mq 



(cf. Swaters et al. 2002 Bajaja et al. 19941 and given the 
total dust mass of 2.9 x 10" M©, the carbon-to-silicate ra- 
tio of 1:3, the LMC gas-phas e C/O ratio of 0.5, the LMC 
gas-phase Si/O ratio of 0.29 fRusse ll fc Dopita|[T992l ) and 
the observed mean oxygen abundance O/H of 1.95 x lO"** 
( [Vigroux et al.|[l987[ ), we estimate that approximately 30 
per cent of carbon, 16 per cent of oxygen and 14 per cent of 
silicon atoms are locked in the solid phase as carbon and sil- 
icate dust. These estimates are consistent with the efflciency 



of 0.12 implied for SN 2003 gd (Sugcrman et al. 20061 and 



the efHciences expected in early-formed galaxies ( [Morgan fc 
Edmundsl|2003 1. They are also broadly consistent with the 



dust condensation efficiency of 0.2 adopted independently 
for the theoretical supernova yields tabulated byWoosley fcl 
Weaver ( [1995) in order to match the dust masses inferred 
for a number of supernovae (Section 3.5.1 1. 



4.3 Spherical symmetry limitations 

To achieve representative heating conditions in a spherical 
model, the numerical setup must allow for enough ionizing 
photons to interact with the correct mass of gas for a given 
strength of the averaged ISRF. The images in Fig. [l] as well 
as the Hubble Space Telescope observations of NGC 4449 



presented by Annibali et al. ( 2008 1 and the study of its 



young stellar clusters by [Reines et aL\ ( |2008[ ) show a varl 
ety of distinct actively star-forming regions. A spherically 
symmetric model must reproduce not only the combined 
strength of the numerous individual sources, but also the 
average ionization structure resulting from the real three- 
dimensional distribution of these ionizing sources. As a con- 
sequence, the required strength of the central ionizing source 
may be slightly overestimated. 

In general. Table |6] shows good agreement between the 
predicted optical and infrared emission line fluxes and the 
observations. Since the ISRF in our models falls with dis- 
tance from the central source as 1/r^, good agreement with 
observations can only be achieved by modifying the dumpi- 
ness or physical distribution of individual atomic or molec- 
ular species. 

For example, our results do not correctly predict the 
ionization structure of sulphur: lines from the higher ion- 
ization potential S^^ ion are systematically underpredicted 
whereas lines from the lower ionization potential ion are 
overpredicted. This discrepancy cannot be explained by the 
assumed global abundances. However, the predicted relative 
sulphur line fluxes shown in Table[6]can be brought into very 
good agreement with observations by using a filling factor 
of £ = 0.055 for sulphur in place of the lower global value 
of 0.033 adopted in this work. Physically, a smoother ISM 
shortens the mean free path of ionizing photons, increasing 
the degree of ionization within the modelled H 11 region. This 
result may suggest that the nebular emission in NGC 4449 
originates from two or more gas phases with differing poros- 
ity, as well as from the transition zones between the phases. 
We conclude that a multi-phase gas model has to be invoked 
in our spherically symmetric approach in order to reproduce 



the details of the observed ionization structure (cf. Cormier 
et al.|2012| ). 



4.4 Distributed ionizing sources 

A more realistic ISRF can be obtained by replacing the sin- 
gle central ionizing source with a uniform distribution of 
sources forming a central 'ionizing sphere'. For this, we as- 
sumed a sphere with r — 0.2 kpc containing 100 identical 
sources, each contributing L*/100 to the total stellar lumi- 
nosity. 

Fig. |6] (top panel, black solid line) presents the predicted 
global SED resulting from the model with distributed ion- 
izing sources. The SED is similar to the best-fitting model 
(red solid line) , except for the range ~ 10-30 fim, where 
the emission is weak in comparison. This discrepancy in the 
strength of the thermal emission by warm dust can be ex- 
plained by the lower average photon energy available to in- 
teract with the dust in the case of distributed sources: dust 
in close proximity to one of the ionizing sources emitting 
with L*/100 is heated to much lower temperatures than is 
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Figure 6. The effects of varying (top panel), the DGR (middle panel) and the relative dust composition (bottom panel) on 
the predicted SED. While L*, the DGR and the distribution of the DGR were free variables in the models, the dust composition 
was fixed (Section |3.5[ ) and its effect is presented here for illustration purposes only. The red solid lines represent the best-fitting 
model. See text for more details. 
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Figure 7. The ionization structure of oxygen as a func- 
tion of radius for the inner 1 kpc in the best-fitting MO- 
CASSIN model of NGC 4449 (dashed lines) and the corre- 
sponding structure for the single ionizing source replaced with 
100 ionizing sources distributed uniformly within a sphere of 
r = 0.2 kpc (solid lines). 

a relatively smaller mass of dust placed near a single source 
emitting with . 

The plots in Fig. [T] illustrate the corresponding differ- 
ences in the ionization structure as a function of radius. 
For example, the overall ionization fraction of O^"*' is signifi- 
cantly higher in the model with distributed sources, whereas 
the overall ionization fraction of O^^ is significantly lower. 
Since in this case the resulting integrated mass of O^"*" is 
smaller, transitions of O''"'' are likely to produce signifi- 
cantly lower line intensities. The opposite would be expected 
for O'^^. Indeed, in the distributed sources model shown in 
Fig. [t] the line fluxes predicted for the [Olll] lines were ap- 
proximately 5 per cent higher, while the flux of the [O ii]a3727 
line was 10 per cent lower than the corresponding fluxes 
in the single source model (Table |6|. As discussed above, 
the key physical factors determining individual line intensi- 
ties are the mean free paths to the ionizing source and the 
strength of that ionizing source. Therefore, in more realis- 
tic models with distributed sources, the input luminosity Li, 
would be expected to be lower by a few per cent to main- 
tain good agreement between the predicted and observed 
line fluxes. 

The physical distribution of ionizing sources affects the 
global ionization structure and, consequently, also affects the 
predicted line intensities, derived elemental abundances (cf. 
Ercolano et al.|2010 1 and may affect the overall stellar lumi- 
nosity L*. Therefore, we note that performing line fitting of 
global emission lines is a valuable constraint on the overall 
properties of a galaxy, but may not be meaningful without 
taking the morphology into consideration. 



5 CONCLUSIONS 

In this paper, we used results from previous studies, the 
available multiwavelength data, the spectral fitting code 
STARLIGHT and a simple chemical evolution code to con- 
struct a photoionization and radiative transfer MOCASSIN 
model of NGC 4449 through an iterative scheme. 

We assumed a simplified star formation history con- 
sisting of three continuous episodes. Using the shape of the 
observed SED, the observed rate of ionizing photons, Q{]1^) 
derived from I/Hq, and the observed nebular emission line 



fiuxes as criteria in our scheme, we infer 3 x 10^ M© for 
the youngest stellar population (0-10 Myr old), which cor- 
responds to an average recent star formation rate of 0.28 
M0 yr~^ or 0.42 M0yr~'^, assuming a Kroupa IMF or a 
Salpeter IMF, respectively. We infer a bolometric stellar lu- 
minosity of 5.7 X 10^ Lq generated by a total stellar mass 
of 1 X 10^ M0. Ahhough the ages of the oldest stellar 
populations in NGC 4449 are not well constrained by the 
model, we note that a very early onset, 12 Gyr ago or ear- 
lier, is selected by an independent analysis using the spectral 
fitting code STARLIGHT. A more in-depth study of stel- 
lar absorption line diagnostics at higher spectral resolutions 
would be necessary to reach more definite conclusions about 
these older populations. 

We modelled the entire stellar component self- 
consistently along with a mixture of carbon and silicate dust. 
Our model yields a dust mass of 2.9 ± 0.5 x 10® Mq, which 
includes 2 per cent of PAHs. This estimate of Mdust could be 
lower, if a higher proportion of carbon-rich dust is assumed. 
Overall, the results from our chemical evolution model for 
a galaxy continuously forming stars are consistent with the 
carbon to silicate dust mass ratio of 1:3 inferred for the LMC 



(Weingartner & Draine 20011. Interestingly, the dust con- 



densation efficiency emerging from our MOCASSIN model 
and from the observed abundance of oxygen in NGC 4449 
is comparable with the value of 0.2 adopted independently 
in our chemical evolution model to match the dust masses 
inferred for a range of supernovae. 

We note that a DGR of 1/190 was derived for the mod- 
elled region with a radius of 3.3 kpc. Including the extended 
H I envelope and the molecular gas is likely to lower the DGR 
to ~ 1/800. Similarly, our model does not account for cold 
dust, which may be responsible for an excess at 850 /im. 
Taking this additional component into account may have a 
significant effect on the derived DGR. 

We conclude that our iterative scheme is a new tool, 
which can be used to model both the dominant stellar pop- 
ulations and the dust content in a self-consistent way. Al- 
though significant degeneracies in the derived parameters 
are expected, they can be reduced by supplementary spec- 
troscopic data. 

In the case of an irregular object, such as NGC 4449, we 
note that the assumption of spherical symmetry may lead 
to a misrepresentation of the ISRF giving rise to unphysical 
conditions near the central ionizing source. In particular, 
our results have shown that a single ISM phase does not 
fully reproduce the observed ionization structure, while the 
assumption of a single DGR overpredicts the mass of warm 
dust, suggesting that the radial density profiles of dust and 
gas are significantly different. 

Our scheme is easily expandable to three dimensions, 
and in future work a realistic distribution of ionizing clusters 
could be used together with a diffuse evolved stellar popula- 
tion component and several gas and dust phases to construct 
more detailed representations of galaxies. In [Paper iT] we 
present a map of the recent star formation rate derived from 
[Nell] 12. 8 + [Ne III] 15. 6, as well as spatially-resolved measure- 
ments of the densities of ionized and neutral ISM phases, 
which may be useful in constructing a three-dimensional 
model of NGC 4449 in the future. 
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APPENDIX A: PHOTOMETRY OF EXTENDED 
SOURCES WITH Swift /XJYOT 



The Ultraviolet/Optical Telescope (UVOT; Roming et al 



2005 1 is one of three instruments on board Swift ( Gehrels 
et al.|[2004[ ) designed to detect and observe gamma-ray 



bursts and their afterglows in seven optical and ultravio- 
let bands. To this end, the instrument is sensitive to single- 
photon events and the data reduction software has been tai- 
lored to point-source observations. However, with its FOV of 
17' X 17', the UVOT can also be of interest for studies of ex- 
tended sources. Although this FOV is considerably smaller 



to that of GALEX (Martin et al. 20051, the UVOT offers 



a higher angular resolution, with PSFs of only 2'.'4-2'.'9 for 
the three UV bands {uvw2, uvm2 and uvwl) covering a com- 
parable wavelength range to one of the two GALEX bands 



(Fig. [Al] [Morrissey et"aL]|2005[ [Breeveld et aL]|2010 l. 

Obtaining reliable photometric measurements for ex- 
tended sources poses a technical challenge, which is inherent 
to the design of the UVOT. The signal from each incoming 
photon is electronically multiplied to generate a splash of 
photons at the CCD stage of the UVOT detector, and the 
centroid of the splash gives positional information to a sub- 
CCD-pixel accuracy. However, the photon-counting detector 
has a readout rate of ~ 90 s~^, which leads to systematic 
undercounting ('coincidence loss' or 'pile-up') when two or 
more photons arrive in a similar location on the detector 
within one frame. In such case, the photons are not only 
counted as one, but the detection is also misplaced by the 
centroiding algorithm. At count rates of ~ 10 s"'^ this effect 



results in a 10 per cent loss in the number of counts ( Poole 



|et al . 2008). 

High-background cases of point-source observations can 
be viewed as close analogues of extended sources. A range of 
background levels was investigated in the models of| Breeveld| 
|et al.| ( |2010[ ), who showed that backgrounds higher than 
~ 0.07 s~ per unbinned image pixel can no longer be fully 
corrected for coincidence loss without introducing an addi- 
tional linear correction factor at each affected pixel. This 
limiting count rate was used by'Hoversten et al. (20111 in 



their analysis of M81 to highlight regions for which coin- 
cidence loss is significant. Those regions in M81 coincided 
with point-like star-forming knots and could be corrected 
individually as point sources in the low background regime. 
However, NGC 4449 is a dwarf starburst galaxy, where the 



star formation activity, comparable to that of M81 ( Gordon 
et al.|2004 \ , is localised in a relatively small volume of space 



The UV count rates in NGC 4449 are high throughout the 
galaxy, as illustrated in Fig. |A2[ making a similar approach 
not feasible. 

Any method of obtaining photometry must take into 
account the variation of emission across the UVOT image, 
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Figure Al. Broad-band effective passbands for UV observa- 
tions with GALEX and Swift. 




Figure A2. Regions affected by coincidence loss in the uvw2- 
band image of NGC 4449. White indicates no coincidence loss. 
North is up, east is to the left. The bar is 1' in length and the 
cross denotes the optical centre of the galaxy. 



as areas with higher count rates suffer from a greater coinci- 
dence loss and require a greater correction. However, coinci- 
dence loss cannot be corrected for on a pixel-by-pixel basis 
because the counts in neighbouring pixels are not indepen- 
dent of each other. The position of each count in the image 
is calculated from a photon splash over five or more physical 
CCD pixels. Consequently, all counts within the same CCD 
pixel (64 image pixels) have been detected through similar 
photon splashes over the same area of the CCD. In the coin- 
cidence loss regime individual detections in a particular area 
on the detector are not independent of each other and should 
be considered collectively for coincidence loss correction. 
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Figure A3. The isophotal setup for the uv'w2-hand image of 
NGC 4449. Four isophotal regions numbered 1—4 are shown in 
grey, red, green and blue, respectively. Region 0, enclosed by 
an ellipse, and a background aperture defined as an elliptical 
annulus, both extend beyond the image and are not shown for 
clarity. The test regions are indicated by 5" circles. North is 
up, east is to the left. The bar is 1' in length and the cross 
denotes the optical centre of the galaxy. 



Al Isophotal correction 

We used NASA's HEAsoft (v6.7; released on 2009-08-19) to 
obtain corrected photometry for three UV and three optical 
bands in NGC 4449. The white band was excluded from this 
study. The task uvotsource was run repeatedly with two 
user-defined apertures: one enclosing a region of interest (or 
'source') and one defining the background. 

The aperture size recommended for UVOT photometry 
is 5" (Poole et al. 20081. By default, if the source aper- 



ture is greater than 5 , the coincidence loss correction fac- 
tor is determined from a 5" circular region at the centre 
of the user-defined aperture. Consequently, for larger aper- 
tures the default procedure is likely to give a significantly 
biased corrected count rate, depending on features present 
in the central area of the user-defined aperture. Therefore, 
it is expected that uvotsource is generally not applicable 
'as-is' to sources of angular extent greater than the size of 
the standard 5" aperture. 

In the method presented below, the image was di- 
vided into regions following suitably chosen isophotal con- 
tours. The coincidence loss correction factors were deter- 
mined from representative 5" 'test' apertures and applied 
to correct photon count rates in the corresponding regions. 
Afterwards, the background was subtracted and the count 
rates were summed to give the total count rate. An illustra- 
tion of this isophotal setup is given in Fig. |A3| Region 1 in 
Fig. I A3| encloses areas with count rates above the threshold 



of 0.028 s per pixel for 2x2 binning, where the effects 



of coincidence loss become non-negligible (Breeveld et al. 
2010[|Hoversten et al.|2011[ ). Regions 2-4 in Fig. |A3| enclose 
areas with count rates of 0.1, 0.2 and 0.3 s~^ per pixel. Fi- 
nally, region is an ellipse enclosing the entire galaxy but 
not aligned with its optical centre. 



Measurements in all individual regions were performed 
by running the HEAsoft command: 

uvotsource iinage=iinage_sk_coadd . f its \ 

expf ile=image_ex_coadd. f its \ 
srcreg=REGION.reg bkgreg=bkg.reg \ 
apercorr=NQNE ceiitroid=NO clobber=ND \ 
frametinie=0 . 0110329 output=ALL sigma=5 \ 
syserr=YES outf il6=photometry .fits 

where REGION. reg is a description of (i) a test region 0-4, 
or (ii) one of the regions comprising regions 0-4 in Fig. |A3[ 

It is useful to define A, Cr(,err), C'test(,err) , /tost and Cbkg(,orr) 

based on the following columns in the output fits table: 



SRC_AREA 

RAW_TOT_RATE(_ERR) 
COI_TOT_RATE(_ERR) 
COI_STD .FACTOR 
COI_BKG_RATE(_ERR) 



area of the source aperture in 
arcsec^, A 

raw count rate (error) in the 
source aperture, Cr(,err) 
corrected count rate (error) in the 
source aperture, Ctest(,err) 
coincidence loss correction factor 
in the source aperture, /tost 

corrected count rate (error) per 
arcsec^ in the background aper- 
ture, Cbkg(,orr). 



The values of Ctcst(,orr) and /tost were measured from 
the five test regions. The measurements from individual re- 
gions in each of regions 0-4 were combined to obtain inte- 
grated measurements for regions 0-4, An, Cr,n and Cr,orr,n. 
Since the same background aperture bkg.reg, defined as 
an elliptical annulus, was used in all instances, Cbkg(,orr) 
represents a single measurement. The total corrected and 
background-subtracted count rate Csrc was then obtained 
by summing individually corrected stripes, or 'isophotes': 



(Al) 



where 



CsTC,n — Ci^n X /tost,ri C'bkg ^ ^i,n, (-^2) 



and 



Ci,, 



Ai^n 



Cr,n — Cr,n + 1 n 3 
Cr 4 n = 4, 



An — An + 1 ^ n ^ 3 

Ai n = 4. 



(A3) 



(A4) 



The upper limit to the statistical error in Csro,n can 
be estimated by assuming that the fractional error in the 
count rate Ci,n in each isophotal region is equal to the frac- 
tional error in the corrected count rates in the corresponding 
test region given by Cteat,err,n/Cteat,n. The errors associated 
with the corrected count rates for higher correction factors 
are significantly higher than the Poisson error, as shown by 

are bino- 



Kuin & Rosen' ( 2008). Both Cto 



and Cbkg 



mial errors associated with the corrected count rates ( Poole 
et al.|2008||Kuin fc Rosen|2008[ ). 
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n 




/tcst,n 




Csrc , n 




[arcsec^] 




[s-1] 


[s-i] 





130000 


1.003 


525.6 ± 1.4 


409 ±21 


1 


8580 


1.034 


470.2 ± 1.1 


479 ±5 


2 


2870 


1.086 


396.8 ±0.9 


428 ±3 


3 


747 


1.134 


181.5±0.5 


205 ± 1 


4 


426 


1.494 


206.8 ±0.2 


309 ± 1 



Table Al. Detailed photometric measurements for the five 
isophotal regions in NGC 4449 in band uvw2 shown in Fig.|A3| 



Cbkg = 9.08 ± 0.03 X 10" 
column definitions. 



* s ^ arcsec ^. See Section 



Al 



for 



This work 



Band 


/max 


C*src 




-^src 


-Fsrc 






[s-i 


] 


[mjy] 


[mJy] 


uvw2/imik 


1.49 


1830 ± 


120 


175 ± 12 


172 ±12 


uvm2 /222lk 


1.32 


1140 ±80 


189 ± 14 


202 ± 14 


uvwl /2486A 


1.43 


1870 ± 


120 


202 ± 15 


207 ±14 


U/3442A 


1.51 


3480 ± 


300 


249 ± 22 


468 ± 40 


6/4321A 


1.40 


4670 ± 


400 


463 ± 41 


1200 ±100 


i;/5410A 


1.18 


2270 ± 


210 


611 ±57 


1250 ±100 



Table A2. Final count rates Csrc and fiuxes Fsrc in six UVOT 
bands for NGC 4449. /m ax gives the maximum correction fac- 
tor used for each band. For comparison, given on the right 
are the fluxes in the 'as-is' approach. All fiuxes have been 
corrected for foreground extinction. 



A2 Results 



Detailed measurements in band uvw2 are presented in T^a.- 
ble |Al[ These results show that even in region 1 the average 
count rate per pixel is 0.055 s~^, which is significa ntly higher 
than the coincidence loss threshold of 0.028 (Hoversten 
et al.|2011| l 



The total count rates in all six UVOT bands are given 
in Table |A2] The corresponding fluxes were obtained using 
the calibrations of |Breeveld| ( |2010| , and were subsequently 
corrected for foreground extinction using E[B — V) = 0.019 
( Schlegcl ct al.^199 8 ) and the extinction law of Cardelli et al 
( 1989| ) with Ry =3.1. 



The uncertainties in the count rates listed in Table [ 
combine (i) the binomial errors in each of the test regions 
0-4 (iKuin I 



Rosen 2008 1 scaled to the area of the corre- 
sponding region, (ii) the estimated uncertainties due to high 
background ( Breeveld et al.|2010 l and (Hi) the uncertainty 
resulting from the choice of thresholds for the five regions. 
For the UV and optical bands, the uncertainties in (ii) were 
estimated at 6 and 8 per cent, respectively. The uncertain- 
ties in (Hi) were found to be at most 2 per cent in all bands. 
Thus, the statistical and systematic uncertainties in the final 
fluxes (Table |A2[ ) amount to ~ 7 and ~ 9 per cent overall 
for the UV and the optical bands. 

Table |A2| also gives total fluxes calculated by UVOT- 
SOURCE directly from region in the 'as-is' approach, i.e., ig- 
noring the spatial variations in the correction factors across 
UVOT images. In bands uvw2, uvm2 and uvwl the fluxes 
agree to within 7 per cent, whereas in the optical bands 
the discrepancy is more significant. The agreement in the 



UV bands may suggest that the centre of region probed 
the global average of the count rate distribution for the UV 
bands, weighted by the extent of coincidence loss. However, 
in sources that are generally fainter in the UV, better agree- 
ment in these bands may be expected as the required correc- 
tion factors are smaller. Although the values of /max listed 
in Table [X2I show that the maximum correction factors were 
high in all bands, the emission in the UV is sharply peaked 
compared to the more uniform emission in the optical. This 
suggests that applying corrections to the UV bands is likely 
to be affected by a smaller systematic uncertainty than ap- 
plying similar corrections to the optical bands. In general, 
these results confirm that uvotsource is not applicable 'as- 
is' to photometric measurements of extended sources, as ex- 
pected from Section [Al| 

The UV and optical fiuxes in Table[2]are consistent with 
those obtained from GALEX and SDSS. A plot of the global 
SED of NGC 4449 in the range 1000-7000 A is presented in 
Fig.jM] 

In the UV, the fluxes derived from GALEX and 
Swift /VYOT agree very well, but the absolute NUV flux is 
slightly lower than the corresponding fluxes obtained from 
Swift /WOT. It should be noted that GALEX, in a similar 
way to Swift /WOT, also suffers from local non-linearities 
near bright sources. These are more difficult to quantify and 
are not automatically corrected for ( Morrissey et al.|2007 |. 
The count rates measured for NGC 4449 in a 3 diameter 
aperture (~ 1000 s"^ in FUV and ~ 4000 s"^ in NUV) sug- 
gest that this effect is likely to lower the measured count 
rates and add a signiflcant systematic uncertainty to the 
derived GALEX fluxes. Therefore, the true GALEX FUV 
and NUV fluxes may be higher and may be associated with 
larger uncertainties than those shown in Fig. |A4[ 

A3 Conclusions 

In the method presented above we used the existing under- 
standing of coincidence loss for point-source observations 
and applied it to obtain photometric measurements for an 
extended source, NGC 4449, in three UV and three optical 
Swift /WOT bands. The derived fluxes are in good agree- 
ment with the available global GALEX and SDSS photom- 
etry, with the overall uncertainties estimated at 7 and 9 
per cent for the UV and the optical bands. 

Extended-source observations with the three narrow 
UV bands and the smaU PSFs of Swift /WOT enable more 
detailed studies of the youngest stellar component in galax- 
ies. Although the uncertainties in our method are compa- 
rable with those of GALEX, we note that at higher count 
rates GALEX also suffers from non-linearity, which is not 
corrected for or included in the published uncertainties. 

This paper has been typeset from a Tp^jX/ M^rjX flle prepared 
by the author. 
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Figure A4. Global photometry of NGC 4449 with GALEX, 
Swift /UVOT and SDSS using data from Table [2] 
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